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Abstract 

The classical Blake threshold indicates the onset of quasistatic evolution leading to cavitation for gas 
bubbles in liquids. When the mean pressure in the liquid is reduced to a value below the vapor pressure, 
the Blake analysis identifies a critical radius which separates quasistatically stable bubbles from those 
which would cavitate. In this work, we analyze the cavitation threshold for radially symmetric bubbles 
whose radii are slightly less than the Blake critical radius, in the presence of time-periodic acoustic 
pressure fields. A distinguished limit equation is derived that predicts the threshold for cavitation for a 
wide range of liquid viscosities and forcing frequencies. This equation also yields frequency-amplitude 
response curves. Moreover, for fixed liquid viscosity, our study identifies the frequency that yields the 
minimal forcing amplitude sufficient to initiate cavitation. Numerical simulations of the full Rayleigh- 
Plesset equation confirm the accuracy of these predictions. Finally, the implications of these findings for 
acoustic pressure fields that consist of two frequencies will be discussed. 

PACS Numbers: Primary 43.25.Yw, Secondary 43.25.Ts, 47.52.-Fj, 43.25.Rq 

Keywords: acoustic cavitation, nonlinear oscillations of gas bubbles, dynamic cavitation threshold, periodic 
pressure fields, quasiperiodic pressure fields, period-doubling. 



I Introduction 



The Blake threshold pressure is the standard measure of static acoustic cavitation j^, |^. Bubbles forced 
at pressures exceeding the Blake threshold grow quasistatically without bound. This criterion is especially 
important for gas bubbles in liquids when surface tension is the dominant effect, such as submicron air 
bubbles in water, where the natural oscillation frequencies are high. 

In contrast, when the acoustic pressure fields are not quasistatic, bubbles generally evolve in highly 
nonlinear fashions ||, |6[ . To begin with, the intrinsic oscillations of spherically symmetric bubbles 
in inviscid incompressible liquids are nonlinear . The phase portrait of the Rayleigh-Plesset equation p6| , 

^, consists of a large region of bounded, stable states centered about the stable equilibrium radius. The 
natural oscillation frequencies of these states depend on the initial bubble radius and its radial momentum, 
and this family of states limits on a state of infinite period, namely a homoclinic orbit in the phase space, 
which acts as a boundary outside of which lie initial conditions corresponding to unstable bubbles. Time- 
dependent acoustic pressure fields then interact nonlinearly with both the periodic orbits and the homoclinic 
orbit. In particular, they can act to break the homoclinic orbit, permitting initially stable bubbles to leave 
the stable region and grow without bound. These interactions have been studied from many points of 
view: experimentally, numerically, and analytically via perturbation theory and techniques from dynamical 
systems. 

In [ p6) , the transition between regular and chaotic oscillations, as well as the onset of rapid radial growth, 
is studied for spherical gas bubbles in time-dependent pressure fields. There, Melnikov theory is applied to 
the periodically- and quasiperiodically-forced Rayleigh-Plesset equation for bubbles containing an isothermal 
gas. One of the principal findings is that, when the acoustic pressure field is quasiperiodic in time with two 
or more frequencies, the transition to chaos and the threshold for rapid growth occur at lower amplitudes of 
the acoustic pressure field than in the case of single-frequency forcing. Their work was motivated in turn by 
that in , where Melnikov theory was used to study the time-dependent shape changes of gas bubbles in 
time-periodic axisymmetric strain fields. 

The work in ||25[| identifies a rich bifurcation superstructure for radial oscillations for bubbles in time- 
periodic acoustic pressure fields. Techniques from perturbation theory and dynamical systems are used to 
analyze resonant subharmonics, period-doubling bifurcation sequences, the disappearance of strange attrac- 
tors, and transient chaos in the Rayleigh-Plesset equation with small-amplitude liquid viscosity and isentropic 
gas. The analysis in ]25t complements the experiments of [^j and the experiments and numerical simulations 
of |l5|, Analyzing subharmonics, these works quantify the impact of increasing the amplitude of the 
acoustic pressure field on the frequency-response curves. 

Other works examining the threshold for acoustic cavitation in time-dependent pressure fields have fo- 
cused on the case of a step change in pressure. In Q, the response of a gas bubble to such a step change in 
pressure is analyzed by numerical and Melnikov perturbation techniques to find a correlation between the 
cavitation pressure and the viscosity of the liquid. One of the principal findings is that the cavitation pres- 
sure scales as the one-fifth power of the liquid viscosity. A general method to compute the critical conditions 
for an instantaneous pressure step is also given in Q. The results extend numerical simulations of |l9| and 
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experimental findings of p4| , and apply for any value of the polytropic gas exponent. 

The goal of the present article is to apply similar perturbation methods and techniques from the theory of 
nonlinear dynamical systems to refine the Blake cavitation threshold for isothermal bubbles whose radii are 
slightly smaller than the critical Blake radius and whose motions are not quasistatic. Specifically, we suppose 
these bubbles are subjected to time-periodic acoustic pressure fields and, by reducing the Rayleigh-Plesset 
equations to a simpler distinguished limit equation, we obtain the dynamic cavitation threshold for these 
subcritical bubbles. 

The paper is organized as follows. In the remainder of this section, the standard Blake cavitation 
threshold is briefly reviewed. This also allows us to identify the critical radius which separates stable and 
unstable bubbles that are in equilibrium. In section ||, the distinguished limit (or normal form) equation of 
motion for subcritical bubbles {i.e., those whose radii is slightly smaller than the critical value) is obtained 
from the Rayleigh-Plesset equation. This necessitates identifying the natural timescale of oscillation of 
such subcritical bubbles which happens to depend upon how close they are to the critical size. We begin 
section [II by defining a simple criterion for determining when cavitation has occurred. We then analyze 
the normal form equation and determine the cavitation threshold for a specific value of the acoustic forcing 
frequency (at which the corresponding linear undamped system would resonate) . This pressure threshold is 
then compared to numerical simulations of the full Rayleigh-Plesset equation and the good agreement found 
between the two is demonstrated. The self-consistency of the distinguished limit equation is further discussed 



in that section. Section IV generalizes the results to include arbitrary acoustic forcing frequencies. Acoustic 
forcing frequencies which facilitate cavitation using the least forcing pressure are determined. An unusual 
dependence of the threshold pressure on forcing frequency is discovered and explained by analyzing the 



"slowly- varying" phase-plane of the dynamical system. At the end of section [V, our choice of a cavitation 
criterion is discussed in the setting of a Melnikov analysis. In section we extend the cavitation results to 
the case of an oscillating subcritical bubble that is driven simultaneously at two different frequencies. We 



recap the paper in section VI by highlighting the main results and discussing their applicability. Lastly, we 
conclude the paper with an appendix which qualitatively discusses the relation of our results to some recent 
experimental findings. 

I.l Blake threshold pressure 

To facilitate the development of subsequent sections we first briefly review the derivation of the Blake 
threshold [0. At equilibrium, the pressure, pb, inside a spherical bubble of radius R is related to the 
pressure, pl, of the outside liquid through the normal stress balance across the surface: 

VB=PL + -^- (1) 

The pressure inside the bubble consists of gas pressure and vapor pressure, pB ~ Pg + Pv, where the vapor 
pressure py is taken to be constant — py depends primarily on the temperature of the liquid — and the 
pressure of the gas is assumed to be given by the equation of state: 

Ps^Psof^) > (2) 
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with 7 the polytropic index of the gas. For isothermal eonditions 7=1, whereas for adiabatic ones, 7 is the 
ratio of constant-pressure to constant-volume heat capacities. At equilibrium, the bubble has radius Rq, the 
gas has pressure Pg^ and the static pressure of the liquid is taken to be Thus, the equilibrium pressure 
of the gas in the bubble is given by 

2cr 

Pgo =Po -P" + 

Upon substituting this result into we get the following expression for the pressure of the gas inside the 
bubble as a function of the bubble radius: 



Upon combining equations (|l|) and (H), we find 



2a\ fRoV'^ 2cr 



Equation (^) governs the change in the radius of a bubble in response to quasistatic changes in the liquid 
pressure pl ■ More precisely, by "quasistatic" we mean that the liquid pressure changes slowly and uniformly 
with incrtial and viscous effects remaining negligible during expansion or contraction of the bubble. For very 
small (sub-micron) bubbles, surface tension is the dominant effect. Furthermore, typical acoustic forcing 
frequencies are much smaller than the resonance frequencies of such tiny bubbles. In this case, the pressure 
in the liquid changes very slowly and uniformly compared to the natural timescale of the bubble. 

For very small bubbles, the Peclet number for heat transfer within the bubble — defined as RqUj/o, with 



Lu the bubble natural frequency (see subsection ll.l) and a the thermal diffusivity of the gas — is small, 
and due to the rapidity of thermal conduction over such small length scales, the bubble may be regarded as 
isothermal. We therefore let 7 = 1 for an isothermal bubble and define 



Then equation (||) becomes 
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The right-hand side of this equation is plotted in figure 1 (solid curve), which shows a minimum value at a 
critical radius labeled -Rcrit- 

Obviously, if the liquid pressure is lowered to a value below the corresponding critical pressure pL„it : 
equilibrium radius exists. For values of which are above the critical value but below the vapor pressure 
equation (^ yields two possible solutions for the radius R. Bubbles whose radii are less than the Blake radius, 
i?crit, are stable to small disturbances, whereas bubbles with R > -Rcrit are unstable to small disturbances. 

The Blake radius itself can be obtained by finding the minimum of the right-hand side of ^ for i? > 0. 
This yields the critical Blake radius 



at which the corresponding critical liquid pressure is 



. (7) 
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By combining the last two equations, it is also possible to express the Blake radius in the form: 



relating the critical bubble radius to the critical pressure in the liquid. Bubbles whose radii are smaller than 
i?crit SLie quasistatically stable, while bigger ones are unstable. 

To obtain the standard Blake pressure we assume that can be ignored and recall that surface tension 
dominates in the quasistatic regime which amounts to ^ 2a/ Rq. Under these approximations, G ~ 2cri?g 
and the Blake threshold pressure is conventionally defined as 

PBlake = - PL„it 

In the quasistatic regime where the Blake threshold is valid, peiakc is the amplitude of the low-frequency 
acoustic pressure beyond which acoustic forcing at higher pressures is sure to cause cavitation. When 
the pressure changes felt by the bubble are no longer quasistatic, a more detailed analysis taking into 
consideration the bubble dynamics and acoustic forcing frequency must be performed to determine the 
cavitation threshold. This is the type of analysis we undertake in this contribution. 



II The distinguished hmit equation 



II. 1 Derivation 



To make progress analytically, we focus our attention on "subcritical" bubbles whose radii are only slightly 
smaller than the Blake radius at a given liquid pressure below the vapor pressure. We thus define a small 
parameter e > by 



e = 2 



1 - 



Rq 



Rc 



(8) 



which measures how close the equilibrium bubble radius Rq is to the critical value i?crit- The value of the 
mean pressure in the liquid, corresponding to the equilibrium radius Rq^ can also be found from equation 
(|) to be 

0{e')]. (9) 
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The liquid pressure and the critical pressure pL„i-t differ only by an 0{e'^) amount. 

It turns out that the characteristic time scale for the natural response of such subcritical bubbles also 
depends on the small parameter e. This timescale for small amplitude oscillations of a spherical bubble is 
obtained by linearizing the isothermal, unforced Rayleigh-Plesset equation 



2cr 



Rq 



2a 

P^-R 



Po 



(10) 



where the density of the liquid is given by p and viscosity has been neglected. Specifically, we substitute 
R = Ro{l + x) into ( p^ ) and keep terms linear in x to get: 



x + 



' Aa 

JrI 



pRl _ 



0. 



(11) 
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Solutions to (p^), representing small amplitude oscillations about equilibrium, are therefore x — xq cos(cjt 
with the angular frequency given by 

_pRl pRl 



(12) 



We now use uj to define a nondimensional time variable: r = ujt. We are interested in analyzing stability for 
values of (i?o,Po°) near (i?crit,PL„it)- Hence, upon recaUing (^, (^ and (^), we see that: 
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crit 



ll/2 



2(Te 



1/2 



We note that as e tends to zero, the timescale for bubble oscillations (the reciprocal of the factor multiplying 
t in the last equation) increases as e~^/^. 

Having determined the proper scaling for the time variable for slightly subcritical bubbles, we can now 
find the distinguished limit (or normal form) equation for such bubbles in a time-periodic pressure field. We 
start with the isothermal, viscous Rayleigh-Plesset equation (2ll: 



RR 



2^ 
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Pv 
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P^+PA sinint). 



(13) 



The amplitude and frequency of the applied acoustic forcing are given by pA and respectively, and 
fi represents the viscosity of the fiuid. Here, the far-field pressure in the liquid has been taken to be 
-Po° ~ Pyi sin(rit), with pg° given by equation (^). Setting R{t) = Ro{l + ex(T)), with e the same small 
parameter introduced above, we obtain at order (noting that all of the 0(1) and 0(e) terms cancel): 



x + 2Cx- 



x^ ^Asin{n*T), 



where 
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(14) 



(15) 



In equation (p^), each overdot represents a derivative with respect to r. 

It is implicit in the above scaling that C, A, and fi* are nondimensional and 0(1) with respect to e. To see 
that this is reasonable, consider an air bubble in water with p — 998 kg/m^, p — 0.001 kg/m-s, a — 0.0725 
N/ni. If we specify e — 0.1 and take a modest equilibrium radius of i?o = 2 x 10^^ m then ( — 0.38. Our 
analysis of (^4|) in subsequent sections will concentrate primarily on values of C, in the range < C < 0.4. 
The parameters A and f2* are related to the forcing conditions, and their magnitudes can be made order 
unity by choosing appropriate forcing parameters pa and il. As an example, if we again choose Rq to 
equal 2 microns, then n* = (2.35 x lO-''s)n/^/e. Moreover, setting e = 0.05 gives n* = (1.05 x lO^^s)!]. 
Hence, the dimensionless parameter fl* is 0(1) when CI is in the megahertz range, and this is precisely 
the frequency range we are interested in exploring. Similarly, with Ro and a chosen as above, we find 
A — (1.38 X lO^^m-s^ /kgJpA/e^ and thereby we see that if e = 0.1 then pA can become on the order of 
lO'^Pa. More data will be presented later, in figure 9, showing typical forcing pressures. 



II. 2 Interpretation 

In the laboratory one can create a subcritical bubble by subjecting the liquid to a low-frequency transducer 
whose effect is to lower the ambient pressure below the vapor pressure. Then a second transducer of high 
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frequency (high relative to the slow transducer) will give rise to the forcing term on the right hand side 
of (^). The low- frequency transducer periodically increases and decreases the pressure in the liquid (and 
shrinks and expands the bubble, which follows this pressure field quasistatically). When the peak negative 
pressure is reached (and the bubble has expanded to its maximum size), we can imagine that state as the 
new equilibrium state, and at that point bring in the effects of the sound from the second transducer. This 
second field can then possibly make the bubble, which had already grown to some large size (but still smaller 
than the critical radius), become unstable. This would all happen very fast compared to the time scale of 
the original slow transducer, so the pressure field contributed by the original transducer remains near its 
most negative value throughout. The stability response of the bubble to the high frequency component of 
the pressure field is the subject of the rest of this work. 



Ill Acoustic forcing thresholds (Q* = 1) 

The value of f7* = 1 corresponds to the forcing frequency at which the linear and undamped counterpart 
of (|lj) would resonate. We therefore choose this value of the forcing frequency as a starting point and 
perform a detailed analysis of the dynamics inherent in the distinguished limit equation at this value of f2*. 
We caution, however, that, as with most forced, damped nonlinear oscillators, the largest resonant response 
occurs away from the resonance frequency of the linear oscillator. We use SI* = 1 mainly as a starting point 



for the analysis, and the dynamics observed for a range of other fl* values is reported in section IV. 

Some special cases of ( p^ ) can be readily analyzed when fl* = 1. In the absence of forcing, i.e., when 
A = 0, the phase portraits of (|l^ ) with C > are shown in figure 2. With no damping (figure 2a), the phase 
plane has a saddle point at (1,0) and a center at (0,0). The latter represents the equilibrium radius of the 
bubble which, when infinitesimally perturbed, results in simple harmonic oscillations of the bubble about 
that equilibrium. The saddle point at (1,0) represents the effects of the second nearby root of the equation 
(I) which is an unstable equilibrium radius. When damping is added (figure 2b), the saddle point remains 
a saddle, but the center at (0,0) becomes a stable spiral, attracting a well-defined region of the phase space 
towards itself. In the presence of weak forcing (small A) but with no damping (C = 0), the behavior of ( p^ ) 
can be seen in a Poincare section shown in figure 3. 



III.l Phase plane criterion for acoustic cavitation 

To determine when a slightly subcritical bubble becomes unstable we choose a simple criterion based upon 
the phase portrait of the distinguished limit equation (|lj). For a given C, there exists a threshold value, 
Acsci of A such that the trajectory through the origin (0,0) grows without bound for A > A^sc, whereas that 
trajectory stays bounded for A < ^osc- A stable subcritical bubble becomes unstable as A increases past 
Acsc- Thus there is a stability curve in the (A, C)-plane separating the regions of this parameter space for 
which the trajectory starting at the origin in the phase-plane either escapes to infinity or remains bounded. 
Numerically, many such threshold C, Acsc pairs (represented by the open circles in figure 4) were found 
with Q* = 1. The data are seen empirically to be well fitted by a least-squares straight line, given by 
Acsc = 1.356C + 0.058. 
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For practical experimental purposes a linear regression curve based upon our escape criterion should 
provide a useful cavitation threshold for the acoustic pressure in the following dimensional form: 

^3/2 1/2 g2^ 

PA > 3.835 ^ + 0.116—. 16 

Here, e is given by equation (^) and is itself a function of the equilibrium radius Rq, surface tension a and 
the pressure differential pg° — pv 

111. 2 Period doubling in the distinguished hmit equation 

It so happens that the stability curve for the trajectory of the origin can also be interpreted in terms of the 
period doubling route to chaos for the escape oscillator (|lj). In other words, the value of A^sc happens to 
be very near the limiting value at which the oscillations become chaotic, just before getting unbounded. For 
a fixed value of C > and a small enough A, the trajectory of the origin will settle upon a stable limit cycle 
in the phase plane. As A is increased gradually, the period of this stable limit cycle undergoes a doubling 
cascade as shown in figure 5 for a fixed value of <^ = 0.35. The period doubling sequence will continue as A 
is increased until the trajectory of the origin eventually becomes chaotic, but still remains bounded. Finally, 
at a threshold value of A the trajectory of the origin will escape to infinity. This is the value of A that is 
given by the open circles on the stability diagram (figure 4). A typical bifurcation diagram for the escape 
oscillator (^^ with C > is shown in figure 6 in which ( = 0.375. 

111. 3 Robustness of the simple cavitation criterion 

In this subsection we justify defining a cavitation criterion based upon the fate of a single initial condition. 
In all simulations with ^ > 0, there is a large region of initial conditions whose fate (escaping or staying 
bounded) is the same as that of the origin (figure 7). In fact, when the trajectory through the origin stays 
bounded it is clear from the simulations that the origin lies in the basin of attraction of a bounded, attracting 
periodic orbit, and points in a large region around it all lie in the basin of attraction of the same orbit. The 
trajectories through all points in that basin remain bounded. Then, as the forcing amplitude is increased the 
attractor is observed to undergo a sequence of period-doubling bifurcations, and this sequence culminates at 
the forcing magnitude when the origin and other initial conditions in a large region about it escape, because 
there is no longer a bounded attracting orbit in whose basin of attraction they lie. 

111. 4 Comparison with full Rayleigh-Plesset simulations 

The stability threshold predicted by equation ( p^ can be compared to data obtained from simulating the 
Rayleigh-Plesset equation ( p^ ) directly. For small values of e, figure 8 shows the resulting good agreement. 

The following is a brief description of how the simulations were carried out. The material parameters 
used to produce figure 8 were p = 998 kg/m"^, /i ~ 0.001 kg/m-s, and cr = 0.0725 N/m. Four values of e were 
chosen, e — 0.01,0.05,0.1, and 0.2. For each fixed value of e and for a selected set of values of ( ranging 
from to 0.4, the parameters Rq, i?crit, [Pv — Po°)j and fl were calculated successively using the formulae: 
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Ro - 2^iyieepa), i?crit = Ro/il-e/2), p.,-p^ - (2a/i?o)[l - (l/3)/(l - e/2)2], and Q = [(2ae)/(pi?3)]i/2. 
(Note that this succession of computations is done for each chosen value of C in each of the four plots.) 

Having obtained the dimensional parameters required for the simulation of the full Rayleigh-Plesset 
equations corresponding to a given (e, Q pair, we used a bisection procedure to determine A^^^ , the threshold 
value of A separating bounded and unbounded bubble trajectories. The bisection procedure was initiated by 
choosing a value of A close to the linear regression line. For this choice of A, the dimensional pressure, pA, 
was calculated using the middle equation in (p^). Then, the initial condition -R(O) = Rq and R{0) — was 
integrated forward in time using an implicit fourth-order Runge-Kutta scheme. The adaptive, implicit 
scheme we used offers an accurate and stable means to integrate the governing equations. The time steps 
are large in those intervals in which the bubble radius does not change rapidly, and they are extremely short 
for the intervals where i? or i? is large (see for example figure 4.7 on page 309 of If the bubble radius 

remained bounded during the simulation, then the value of A was increased slightly, a new pA was calculated, 
and a new simulation was begun. If, on the other hand, the bubble radius became unbounded during the 
simulation, then the value of A was slightly decreased and a new simulation was initiated. Continuing with 
this bisection of A, the threshold value, A^^^ , where the bubble first becomes unstable was determined. 

The dimensional counterpart (pA-'versns-Ro) to figure 8 is shown in figure 9 along with the dimensional 
stability curve given by equation (p^. Note that for a given parameter e, the relationship which defines the 
dimensionless damping parameter C, i.e., Rq = 2/i'^/(C^epcr), can be thought of as defining the bubble radius 
_Ro- That is, for a given liquid viscosity /i and with all other physical parameters being constant, C can only 
change by varying the equilibrium radius Rq. As such, the dimensionless A-versus-C curves can be put in 
terms of the dimensional pA-veisus-Ro curves, drawn in figure 9. 

To show the way in which the bubble radius actually becomes unbounded in the full Rayleigh-Plesset 
simulations, figure 10 provides the radius- versus-time plots for three typical simulations with the same value 
of C = 0.3, where time is nondimensional. In this case, A^^^ w 0.51. The top two curves are obtained 
for values of A of 0.3 and 0.5, respectively. They show stable oscillations although a period-doubling can 
be seen to have occurred in going from one to the other. The bottom figure corresponds to A = 0.53 and 
shows that the bubble radius is becoming unbounded. The corresponding dimensional parameters for the 
Rayleigh-Plesset simulations are given in the figure caption. 

III. 5 Consistency of the distinguished hmit equation 

In this subsection, we argue, a posteriori, that it is self-consistent to use the escape oscillator as the dis- 
tinguished limit equation ( [T^ ) for the full Rayleigh-Plesset (^3|), i.e., we show that the higher-order terms 
encountered during the change of variables from R to x may be neglected in a consistent fashion. 

Recall that, in the derivation of the escape oscillator, all of the 0(1) and 0{e) terms dropped out, and 
the ordinary differential equation ( [T^ was obtained by equating the terms of 0{e^). The remaining terms 
are of ©(e^) and higher. To be precise, at C(e^), we find on the left-hand side: 

XX ^(^XX I X ^ 
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and on the right-hand side: 

3 „ , 20 , 

4 3 

Moreover, we note that, for i > 4, aU terms of 0{e^) on the left-hand side are of the form a;*~^i;, while all 
terms of O(e') on the right-hand side are polynomials in x. 

We already know that, for trajectories of the escape oscillator that remain bounded, the x and x vari- 
ables stay 0(1). Hence, all of the higher-order terms remain higher-order for these trajectories. Next, for 
trajectories that eventually escape (i.e., those whose x-coordinate exceeds some large cut-off at some finite 
time), we know that x and x are bounded until that time and afterwards they grow without bound. The 
fact that then R also grows without bound for these trajectories (due to the change of variables that defines 
x) is consistent with the dynamics of the full Rayleigh-Plesset equation. 

Potential trouble could arise with trajectories for which x becomes negative and large in magnitude, 
e.g., when x ^ — 1/e the coefficient of x vanishes. (This corresponds to small R.) A glance at the Poincare 
map for the escape oscillator reveals, however, that trajectories which have x ~ — 1/e at some time r can 
never have x{t) and x{t) of 0(1) simultaneously, for any t. Hence, these trajectories are not in the regime 
of interest, neither for the escape oscillator nor for the full Rayleigh-Plesset equation. This completes the 
argument that it is self-consistent to use the escape oscillator for this study. 



IV Pressure thresholds for general Q* 
IV. 1 Stability curves for various Q* 

Until now, we have only examined the case il* = 1 in the distinguished limit equation (|l^). In this section 
we examine the dependence of the stability threshold on the acoustic frequency, fl, for subcritical bubbles. 
Specifically, for frequencies Q* between 0.1 and 1.1, we performed numerical simulations of the distinguished 
limit equation ( p^ to determine many (C, ^osc) pairs. These pairs are plotted in figure 11, and the data 
points at each dimensionless frequency fl* are connected by straight lines (in contrast to the least squares 



fitting done in subsection [II.l). 

As in subsection III.l, good agreement between the distinguished limit equation threshold and the full 
Rayleigh-Plesset equation is observed for various values of fi*; this can been seen in figure 12 which compares 
the two results at four different values of H.* given by 0.6, 0.7, 0.8 and 0.9, for a fixed value of e = 0.05. 

Various features observed in figure 11, such as the flattening of these curves as ft* decreases, are explained 
in the next subsection. 



IV. 2 Minimum forcing threshold 

Suppose that we wish to determine the driving frequency, for a given bubble with an equilibrium radius _Ro 
and a critical radius i?crit, so that the acoustic forcing amplitude necessary to make the bubble unstable is 
minimized. This can be done by choosing in figure 11, the value of Q* for which the corresponding threshold 
curve is below all the others for a given ^. The result of such a procedure is provided in figure 13 as follows. 
Figure 13(a) provides the frequency of harmonic forcing at a given value of the damping parameter C for 
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which the required amphtude of the acoustic pressure field to create cavitation is the smallest. Figure 13(b) 
shows the dimensionless minimum pressure amplitude ^csc corresponding to the value of il* just presented. 
In figure 13(a), for ^ between and 0.225, the frequency curve is nearly a straight line and we fit a linear 
regression line to the data in that interval: fi* = —1.12^ + 0.90 for < ^ < 0.225. Correspondingly, in figure 
13(b), we see that the minimum pressure curve is also nearly straight for the same interval of ^ values. The 
least squares line fitting the data in figure 13(b) is ^ = 1.03^ + 0.02 for < < 0.225. 

When ( — 0.225, there is a discontinuity in the frequency curve, as seen in figure 13(a). At the same 
value of C, the pressure curve levels off to yl sa 0.25. This can be explained by a brief analysis of the normal 
form equation. The key observation will be that, in the escape oscillator with constant forcing (i.e., constant 
right-hand side), there is a saddle-node bifurcation when the magnitude of the forcing is 1/4. In order to 
carry out this brief analysis, we consider the cases C > ^-i^d C = separately, beginning with C > 0. 

For C > 0, the Poincarc map of the normal form equation (|l^ ) has an asymptotically stable fixed point 
(a sink), which corresponds to an attracting periodic orbit for the full normal form equation. Now, during 
each period of the external forcing, the location of this periodic orbit in the (x, i;)-plane changes. In fact, 
for the small values of Q* we are interested in here (f2* < 0.3 approximately), the change in location occurs 
slowly, and one can write down a perturbation expansion for its position in powers of the small parameter 
51*. The coefficients at each order are functions of the slow time z = il*t. To leading order, i.e., at 0(1), the 
attracting periodic orbit is located at the point (a;(z),0), where x{z) is the smaller root oix — x^ — Asin(z), 
namely 

x{z) = i - iv/l-4Asin(z). 

Therefore, one sees directly that A = 1/4 is a critical value. In particular, if one considers any fixed 
value of ^ < 1/4, then the attracting periodic orbit exists for all time, and the trajectory of our initial 
condition (0, 0) will be always be attracted to it. (Note that the viscosity Q > 0.225 is large enough 
so orbits are attracted to the stable periodic orbit at a faster rate than the rate at which the periodic 
orbit's position moves in the {x, a;)-plane due to the slow modulation.) However, for any fixed value of 
A > 1/4, the function giving x{z) becomes complex after the slow time z reaches a critical value 2* (A), 
where Asm{z^) ~ 1/4, and where we write z^{A) since depends on A. Moreover, x{z) remains complex 
in the interval (z*(yl), tt — z*(A)) during which j4sin(z) > 1/4. Viewed in terms of the slowly-varying phase 
portrait, the slowly- moving sink merges with the slowly-moving saddle in a saddle-node bifurcation when z 
reaches z^{A), and they disappear together for z*(A) < z < n — z^{A). Hence, the attracting periodic orbit 
no longer exists when z reaches z^{A), and the trajectory that started at (0,0) — and that was spiraling in 
toward the slowly-moving attracting periodic orbit while z was less than z^,{A) — escapes, because there is 
no longer any attractor to which it is drawn. 

For the sake of completeness in presenting this analysis, we note that when A = 1/4, then z^{A) — 7r/2; 
hence, it is precisely near this lowest value of A, namely A — 1/4, that we find the threshold for the acoustic 
forcing amplitude, and the escape happens near the slow time z — 7r/2. Moreover, for values oi A > 1/4, 
< z^{A) < tt/2, and so the escape happens at an earlier time. 

Numerically, the minimal frequency fl* appears to be H.* « 0, where ft* = 0.01 is the lowest value for 
which we conducted simulations. Moreover, this also explains why, as we see from figure 11 already, the 
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curves are fiat with A « 0.25 for H.* < 0.3. Tliis is tlie range of small values of fl* for which the above 
analysis applies. 

Next, having analyzed the regime in which C > 0, wc turn briefly to the case ( = 0- For small values, 
fl* < 0.3, the curves in figure 11 remain flat near A = 0.25 all the way down to = 0. The full normal form 
equation with C — is a slowly-modulated Hamiltonian system. One can again use the slowly-varying phase 
planes as a guide to the analysis (although the periodic orbit is only neutrally stable when C = and no 
longer attracting as above), and the saddle- node bifurcation in the leading order problem at Asm{z^) = 1/4 
is the main phenomenon responsible for the observation that the threshold forcing amplitude is near 0.25. 
(We also note that for a detailed analysis of the trapped orbits, one needs adiabatic separatrix-crossing 
theory, see for example, but we shall not need that here.) 

Finally, and most importantly, simulations of the full Rayleigh-Plesset equation confirm all the quan- 
titative features of this analysis of the normal form equation. The open circles in figure 13 represent the 
numerically observed threshold forcing amplitudes, and these circles lie very close to the curves obtained as 
predictions from the normal form equation. We attribute this similarity to the fact that the phase portrait 
of the isothermal Rayleigh-Plesset equation has the same structure — stable and unstable equilibria, sep- 
aratrix bounding the stable oscillations, and a saddle-node bifurcation when the forcing amplitude exceeds 
the threshold — as the normal form equation (see p^). 

IV. 3 The dimensional form of the minimum forcing threshold 

Recall, that for C between and 0.225, we fit linear regression lines to portions of figures 13(a) and (b). 
Specifically, for figure 13(a) we found that, for a particular choice of (, the frequency which yields the smallest 
value of Aesc can be expressed as: il* = —1. 12^ + 0. 90 for < C < 0.225. And for figure 13(b) we found that 
the stability boundary for the minimum forcing is given by. 



Using the definitions of (,A and Q* as given by (|l5|), the "optimal" acoustic frequency to cause cavitation 
of a subcritical bubble is given in dimensional form by 




for < C < 0.225 
for 0.225 <C<0.4. 



(17) 




( 



2^2 



) 



1/2 



n = -2.24-^ + 1.27 



for 



< 



< 0.225. 



aepRo 



Correspondingly, the minimum acoustic pressure threshold is 




IV. 4 A lower bound for ^csc via Melnikov analysis 



The distinguished limit equation dT3) can be written as the perturbed system 



x = f(x) +eg(x,r) 
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where, x — {x,y), i{x,y) — [y^x^ — x) and ^[x,y,T) = (O, Asin(r2*T) — 2(^1/) with A = lA, ( — e(. 
When e = the system has a center at (0,0) and a saddle point at (1,0). The homocHnic orbit to the 
unperturbed saddle is given by 7o(t) = {x{T),y{T)) where x{t) = -(1/2) + (3/2) tanh^(T/2) and ^(t) = 
(3/2)tanh(r/2)sech^(T/2). Following jll, the Melnikov function takes the form 



M(to) = / f(7o(T)) Ag(7o(T),T + To)dT 

J — oo 

= ^A f sin[f7*(T + To)] tanh sech^ dr 

-^C/Ja.l,^(0»ch«(l)... 

The first integral can be done with a residue calculation, and the second integral is evaluated in a straight- 
forward manner, resulting in: 



"64f]*)2cos(r!Vo) 



sinh(7rri*) 

The Melnikov function has simple zeros when A > ^h.tan., where 



_ , 2sinh(7rl]*) , 

^h.tan. - ( 5^(^,)2 ) C- (18) 



Hence the stable and unstable manifolds of the perturbed saddle point intersect transversely for all sufficiently 
small e 7^ when A > Ah. tan. l|l2|. The resulting chaotic dynamics is evident in figure 5, for example. Since 
homoclinic tangency must occur before the trajectory through the origin can escape, Ah. tan. may be viewed 
as a precursor to A^sc- Figure 14 demonstrates that, for small enough e, equation ( p^ ) provides a lower 
bound for the stability curves seen in figure 11. 

The reason why Melnikov analysis yields a lower bound for the cavitation threshold relates to how 
deeply the stable and unstable manifolds of the saddle fixed point of the Poincare map for equation ( p^ ) 
penetrate into the region bounded by the separatrix in the A, C = case. For sufficiently small values of e, 
long segments of the perturbed local stable and unstable manifolds will stay 0{e) close to the unperturbed 
homoclinic orbit. However, as e grows (and one gets out of the regime in which the asymptotic Melnikov 
theory strictly applies), these local manifolds will penetrate more deeply into the region bounded by the 
separatrix in the A,(^ — case. In fact, there is a sizable gap in the parameter space between the homoclinic 
tangency values and the escape values corresponding to our cavitation criteria, i.e., when the trajectory 
through the origin grows without bound. There is a similar gap when other initial conditions are chosen. 

The Melnikov function was also calculated in There, a detailed analysis of escape from a cubic 

potential is described and the fractal basin boundaries and occurrence of homoclinic tangencies are given. 
We also refer the reader to jl^ in which a closely related second order, damped and driven oscillator with 
quadratic nonlinearity is studied using both homoclinic Melnikov theory, as was done here, and subharmonic 
Melnikov theory. The existence of periodic orbits is demonstrated there, and period doubling bifurcations 
of these periodic orbits are examined. Their equation arises from the study of travelling waves in a forced, 
damped KdV equation. 
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V Pressure fields with two fast frequencies 



In this section, we consider what happens to the cavitation threshold if two fast frequency components 
are present in the acoustic pressure field, and the slow transducer, which lowers the ambient pressure and 
whose effect is quasistatic, is also still present. In figure 15, we show the results from simulations with 
quasiperiodic pressure fields. These were obtained from simulations of ( p^ ) with the forcing replaced by 
(A/2)(sin(rijT) + sin(ri2T)), and a wide range of values for fll and il^. For a fixed value of C = 0.25, the 
cavitation surface shown in the figure was plotted by computing the triples (^11, i^2i ^csc)- 

We note that in figure 15, the intersection of the cavitation surface and the vertical plane given by 
ill — ^2 represents cavitation thresholds for acoustic forcing of the form Asm{il*T) {i.e., a single fast 
frequency component and a quasistatic component). Furthermore, we see that the global minimum of the 
cavitation surface lies along the line fJ* — ilj. Hence, for A/2 as our particular choice of quasiperiodic 
forcing coefficient, the addition of a second fast frequency component in the pressure field does not lower 
the cavitation threshold beyond that of the single fast frequency case. 

VI Discussion 

A distinguished limit equation has been derived which is suitable for use in determining cavitation events of 
slightly subcritical bubbles. This "normal form" equation allows us to study cavitation thresholds for a range 
of acoustic forcing frequencies. For fi* = 1, we find an explicit expression for the cavitation threshold via 
linear regression, since the simulation data reveal an approximate linear dependence of the nondimensional 
threshold amplitude. A, on the nondimensional liquid viscosity, C- When converted to dimensional form, this 
linear expression translates into a nonlinear dependence, cf. equation (^|), on the material parameters. In 
all of our simulations, the acoustic threshold amplitude coincides with the amplitudes at which the cascades 
of period-doubling subharmonics terminate. 

Particular attention has also been paid to calculating the frequency, fl* , at which a given subcritical 
bubble will most easily cavitate. Expression (^|) for the corresponding minimum threshold amplitude A 
grows linearly in ^ for ^ < 0.225 until the critical amplitude A = 1/4 is reached, and the threshold amplitude 
stays constant at A w 1/4 for larger C- For these larger values of C > 0.225, the "optimal" frequency is 
essentially zero, as we showed by doing a slowly- varying phase portrait analysis and exploiting the fact that 
the normal form equation undergoes a saddle- node bifurcation at A ^ 1/4 in which the entire region of 
bounded stable orbits vanishes. The full Rayleigh-Plcsset equation undergoes a similar bifurcation at forcing 
amplitudes very near ^ = 1/4 for sufficiently small e. Overall, the results from the normal form equation 
are in excellent agreement with those of the full Rayleigh-Plesset equation, and this may be attributed to 
the high level of similarity between the phase-space structures of both equations. 

In view of the findings in [ p6| , we may draw an additional conclusion from the present work. In a certain 
sense, we have extended the finding of lowered transition amplitudes reported in p6| to the limiting case 
of one low frequency and one fast frequency. We find that if a low frequency transducer prepares a bubble 
to become slightly subcritical, then the presence of a high frequency transducer can lower the cavitation 
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threshold of the bubble below the Blake threshold. 

Our results on the optimum forcing frequency and minimum pressure threshold to cavitate a subcritical 
bubble may also be useful in fine-tuning experimental work on single-bubble sonoluminescence (SBSL). In 
SBSL jl^, ^, |2^, a single bubble is acoustically forced to undergo repeated cavitation/coUapse cycles, in 
each of which a short-lived flash of light is produced. While the process through which a collapsing bubble 
emits light is very complex and involves many nonlinear phenomena, the possibility of better control over 
cavitation and collapse, e.g., through the use of multiple-frequency forcing, can perhaps be investigated using 
the type of analysis presented in this paper. 
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Appendix: Coaxing experiments in acoustic microcavitation 

To illustrate one application of our results, we now briefly consider the experimental findings of p8[ | on 
so-called "coaxing" of acoustic microcavitation. In these experiments, smooth submicrometer spheres were 
added to clean water and were found to facilitate the nucleation of cavitation events (i.e. reduce the cavitation 
threshold) when a high-frequency transducer (originally aimed as a detector) was turned on at a relatively 
low pressure amplitude. Speciflcally, the main cavitation transducer was operating at a frequency of 0.75 
MHz, while the active detector had a frequency of 30 MHz. In a typical experiment, with 0.984-micron 
spheres added to clean water, the cavitation threshold in the absence of the active detector was found to be 
about 15 bar peak negative. When the active detector was turned on, producing a minimum pressure of only 
0.5 bar peak negative by itself, it caused the cavitation threshold of the main transducer to be reduced from 
15 to 7 bar peak negative. The polystyrene latex spheres were observed under scanning electron microscopes 
and their surface was determined to be smooth to about 50 nanometers. It was thus thought that any gas 
pockets which were trapped on their surface due to incomplete wetting and which served as nucleation sites 



for cavitation, were smaller in size than this length. In |18 , it is conjectured that the extremely high fluid 
accelerations created by the high-frequency active detector, coupled with the density mismatch between the 
gas and the liquid, caused these gas pockets to accumulate on the surface of the spheres and form much 
larger "gas caps" (on the order of the particle size), which then cavitated at the lower threshold. Here we 
attempt to provide an alternative explanation for the observed lowering of the threshold in the presence of 
the active detector. 

To effect our estimates, we shall use the same physical parameters as earlier: — 0.001 kg/m-s, p ~ 0.998 
kg/m^^ and <j ~ 0.0725 N/m. We also ignore the vapor pressure of the liquid at room temperature. Note also 
that 1 bar=10^ N/m'^ and that the transducer frequencies / cited above are related to the radian frequencies 
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fl used earlier hy fl = 27r/. 

Let us begin by estimating a typical size for the nucleation sites which cavitate at pLcit = bar in the 
absence of the active detector. Upon using Blake's classical estimate oi pl„h — — 0.77o'/i?o, the equilibrium 
radius of the trapped air pockets is estimated to be Rq — 3.7 x 10^® m or 37 nm. This size is consistent with 
the observation that the surfaces of the spheres were smooth to within 50 nm. We note that such a small 
cavitation nucleus cannot exist within the homogeneous liquid itself since it would dissolve away extremely 
fast due to its overpressure resulting from surface tension. However, when trapped in a crevice or within 
the roughness on solid surfaces, it can be stabilized against dissolution with the aid of the meniscus shape 
which separates it from the liquid. The natural frequency of a 37 nm bubble (if it were spherical) found 
from equation (|2|) would be 385 MHz which is very large compared to forcing frequency of the cavitation 
transducer which is 0.75 MHz. Therefore, consistent with Blake's classical criterion, the pressure changes in 
the liquid would appear quasistatic to the bubble and at such a small size, surface tension does dominate 
the bubble dynamics. The Blake critical radius i?crit which corresponds to this equilibrium radius Rq of 37 
nm can be calculated to be i?crit — 64 nm. 

Let us now suppose that the cavitation transducer is operating at pl = —7 bar peak negative as in the 
experiments with the active transducer also turned on. Using equation , the final expanded radius of the 
bubble when the liquid pressure is quasistatically reduced to —7 bar is found to be 4.2 x 10^^ m or 42 nm. 
In other words, a bubble of original radius 37 nm at a liquid pressure of 1 bar, grows to a maximum size of 
42 nm when the liquid pressure is reduced to —7 bar. Its critical radius is still 64 nm, reached if the liquid 
pressure were to be reduced further to —15 bar. 

At this point, since the pressure changes in the liquid due to the 0.75 MHz cavitation transducer are 
occurring slowly compared both with the natural timescale of the bubble and the 30 MHz detector, let us 
take the mean pressure in the liquid to be the ~ —7 bar, and imagine the bubble size at this pressure 
to be its new equilibrium radius Rq = 42 nm, with the critical radius still given by i?crit = 64 nm. This 
bubble is now assumed to be forced by the 30 MHz transducer at an acoustic pressure amplitude of 1.5 bar 
{i.e. —0.5 bar peak negative). Using these values, the perturbation parameter e is calculated from equation 

to be e = 0.69. This parameter is too big for the results of the asymptotic theory to provide meaningful 
quantitative agreement; nevertheless, we proceed with the discussion to see if we can at least obtain the 
right order of magnitude for the pressure threshold. 

With the given physical parameters, and using the forcing pressure of pa = 1.5 bar and fl — 27r x 30 x 10~^ 
s^^, the parameters C, ^ and fl* are calculated from equation ( |l5| ) to be: ( = 0.98, A = 0.09 and fl* ~ 0.16. 
Upon examining figure 13, at the relatively large damping parameter C = 0.98 (beyond the range originally 
considered) the minimum forcing threshold would appear to correspond to the constant value oi A = 0.25. 
Here we also note that this same forcing threshold is also observed with a range of small fJ*, including 
ri* = 0.16, see figure 11. In other words, the predicted threshold pressure for the active detector to cause 
cavitation is A = 0.25 which corresponds roughly top a =4 bar, whereas in the experiments the threshold was 
seen to be A = 0.09 or = 1.5 bar. Despite the lack of quantitative agreement, the theoretical predictions 
and the experiments do show the same trends. Namely, in the absence of the 30-MHz detector, the pressure 
in the liquid had to be reduced to —15 bar for cavitation to occur. With the high-frequency transducer 
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turned on, however, cavitation occurred at a minimum pressure of —7 — 1.5 = —8.5 bar in the experiments 
and at —7 — 4= —11 bar based on the theory. (We are adding the negative pressure contribution from the 
two transducers to arrive at the final minimum pressure). Thus, the presence of the second high-frequency 
transducer does reduce the pressure threshold for cavitation in both cases. 
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Figure Captions 



Figure 1: Pressure in the liquid, p^, versus bubble radius, R, as governed by equation (|5|). 

Figure 2: Phase portraits for the distinguished limit equation (|lj). In (a), ^ = and C = 0. The fixed 
point (0,0) is a center. The fixed point (1,0) is a saddle. In (b), A = and C — 0.09. The fixed point (0,0) 
is a stable spiral. 

Figure 3: Poincare section showing the unstable manifold of the saddle fixed point (0.999769375,-0.024007197) 
for A — 0.048 and C = 0. Asymptotically, the saddle point is located at a distance 0{A) from (1,0), specif- 
ically (1 - A^/IO + 0{A^), -A/2 - (13/200)^3 -|- ©(A^)). Invariant tori are shown inside a portion of the 
unstable manifold. The center point has moved a large distance from (0,0) in this case due to the 1:1 
resonance when 51* = 1. And we note for comparison that with nonresonant values of fi*, Poincare sections 
show that the center only moves an 0{A) distance. For example, when fl* = 0.6,0.7,0.85, the centers are 
located approximately at the points (0.009,0.045), (0.009,0.066) and (0.029,0.163) respectively. 

Figure 4: Escape parameters for the trajectory of the origin (0,0). For values of A,( above the regression 
line the trajectory of the origin grows without bound. Below this line the trajectory of the origin remains 
bounded. Least squares fit: A = l.356( + 0.058. 

Figure 5: Period doubling route to chaos in the distinguished limit equation (|lj). For C — 0.35, the limit 
cycles undergo period doubling as A is increased. 

Figure 6: Bifurcation diagram {( = 0.375). Plotted is i versus A. For each fixed value of A, the origin 
(0,0) is integrated numerically and the value of x is plotted every At — 27r. 

Figure 7: Bounded trajectories for the distinguished limit equation with ( — 0.2, A — 0.3 and il* = 1.0. 
The dark region is the set of initial conditions whose trajectories remain bounded; it is the basin of attraction 
of the periodic orbit that exists in the period-doubling hierarchy for this value of A. 

Figure 8: Simulations of the full Rayleigh-Plesset equation for four different values of e. Each open 
circle represents an {A, C) pair at which the bubble first goes unstable. Superimposed is the linear re- 
gression line obtained from the simple criterion based upon the distinguished limit equation. In (a)-(d), 
e = 0.01, 0.05, 0.1, 0.2, respectively. 



Figure 9: Simulations of the full Rayleigh-Plesset equation for four different values of e. Each open 
circle represents a (pA, Ro) pair at which the bubble first goes unstable. Superimposed is the threshold 
curve (|l6|) obtained from the simple criterion based upon the distinguished limit equation. In (a)-(d), 
e = 0.01,0.05,0.1,0.2, respectively. 



Figure 10: Radius versus time plots (e = 0.1, C = 0.3, A = 0.3, 0.5, 0.53). 

Dimensional parameters: Rq ~ 3.0 /im, Rcrit = 3.2 /im, p„ — = 29.7 kPa, fl — 0.7 MHz. 

Top: Pa = 141.6 Pa. Middle: pa = 236.0 Pa. Bottom: pA = 250.2 Pa. 



Figure 11: Stability threshold curves for the origin trajectory of the distinguished limit equation for 
many different values of fi* . The values of f2* on the right hand border label the different threshold curves. 



Figure 12: Simulations of the full Rayleigh-Plesset equation for four different values of fl*. In each plot 
e — 0.05. Each open circle represents an A, ^ pair at which the bubble first goes unstable. Superimposed 
is the threshold curve obtained from the distinguished limit equation. In (a)-(d), fi* — 0.6,0.7,0.8,0.9, 
respectively. 



Figure 13: In (a), the il* that minimizes Acsc is plotted versus (. In (b), the minimum value of ylcsc 
corresponding to the value of fl* in (a) is plotted versus (■ The open circles represent Rayleigh-Plesset 
calculations of the minimum threshold with e = 0.05. 



Figure 14: Comparisons of the stability curves of figure 11 with the Melnikov analysis for two values of fl* . 
The dotted line is the stability curve obtained from the distinguished limit equation. The solid straight line 



is equation (18). In (a) and (b), fi* = 0.9, 1.1, respectively. 



Figure 15: Stability threshold surface for the trajectory through the origin obtained by integrating the 
distinguished limit equation with quasiperiodic forcing. This was obtained from simulations of (|lj) with 
the forcing replaced by {A/2){sm{illT) + sin(f72T)), and for < ill,n2 < 1.3. The value of ( was fixed at 
C = 0.25. The points below the surface correspond to parameter values for which the trajectory of the origin 
remains bounded whereas those points above the surface are parameters which lead to an escape trajectory 
for the origin. Qualitatively and quantitatively similar results were obtained for ^ = 0.15 and ^ — 0.35. 
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